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Abstract 

Preconditioned iterative methods for the indefinite systems obtained by dis- 
cretizing the linear elasticity and Stokes problems with mixed spectral elements 
in three dimensions are introduced and analyzed. The resulting stiffness matrices 
have the structure of saddle point problems with a penalty term, which is associ- 
ated with the Poisson ratio for elasticity problems or with stabilization techniques 
for Stokes problems. The main results of this paper show that the convergence 
rate of the resulting algorithms is independent of the penalty parameter, the num- 
ber of spectral elements N and mildly dependent on the spectral degree n via the 
inf-sup constant. The preconditioners proposed for the whole indefinite system are 
block-diagonal and block-triangular. Numerical experiments presented in the final 
section show that these algorithms are a practical and efficient strategy for the 
iterative solution of the indefinite problems arising from mixed spectral element 
discretizations of elliptic systems. 


*This work was supported by the National Science Foundation under Grant NSF-CCR-9503408 and 
by the National Aeronautics and Space Administration under NASA Contract No. NASl-19480 while 
the author was in residence at the Institute for Computer Applications in Science and Engineering 
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1 Introduction 


The goal of this paper is to introduce and analyze some preconditioned iterative meth- 
ods for the large indefinite linear systems arising from the mixed spectral discretization 
of the linear elasticity system and the limiting Stokes problem in three dimensions. 
Standard finite element discretizations of elasticity problems can suffer from the phe- 
nomenon of locking when the Poisson ratio tends to 1/2 (almost incompressible case); 
see Babuska and Suri [3]. This means that the convergence rate of the finite element 
method deteriorates when v approaches 1/2. Moreover, the resulting linear system, 
even though symmetric and positive definite, has a condition number that goes to in- 
finity when the Poisson ratio tends to 1/2. Both problems can be overcome by using a 
mixed finite element formulation and rewriting the problem as a saddle point problem 
with a penalty term; see Brezzi and Fortin [12]. The penalty term depends on the Pois- 
son ratio for elasticity problems or on stabilization parameters for Stokes problems. By 
carefully choosing the finite element spaces in order to satisfy the inf-sup condition, we 
obtain a convergent method. The stiffness matrix is symmetric and indefinite. 

In recent years, several iterative methods have been proposed and studied in the 
case of low-order h-version finite elements, such as Uzawa’s algorithm (see Elman and 
Golub [17], Bramble, Pasciak and Vassilev [10]), multigrid (see Verfiirth [37], Wit- 
tum [39], Braess and Blomer [8], Brenner [11]), preconditioned conjugate gradient (see 
Bramble and Pasciak [9]), and preconditioned conjugate residuals (PCR) (see Rusten 
and Winther [34], Silvester and Wathen [35], [38], Klawonn [22], [24]). Elman [16] has 
carried out a careful comparison of the performance of four of these methods applied 
to Stokes problems in two dimensions. 

Here we consider instead spectral element discretizations. For a general introduc- 
tion to spectral methods, we refer to Canuto, Hussaini, Quarteroni, and Zang [13], 
Bemardi and Maday [7], and Funaro [20]. See also Babuska and Suri [4] for the re- 
lated p - version of the finite element method. Already for scalar problems, the stiffness 
matrices obtained by spectral and /> version finite elements are less sparse and more 
ill-conditioned than those obtained with /i- version finite elements. The construction 
and analysis of efficient preconditioned iterative methods is therefore more challeng- 
ing. We refer to Pavarino and Widlund [31], Casarin [14] and to the references therein 
for an overview of recent results based on domain decomposition techniques for elliptic 
scalar problems. In the context of spectral elements for Stokes and Navier- Stokes prob- 
lems, iterative methods have been studied in Maday, Meiron, Patera and Rpnquist [25], 
Maday, Patera and Rpnquist [26], Fischer and Rpnquist [18], and Rpnquist [33]. The 
methods proposed by these authors are based oh conjugate gradient iterations on the 
reduced Schur complement of the discrete Stokes matrix involving only the pressure 
unknowns. In the context of linear elasticity and p- version finite elements, iterative 
methods have been studied by Mandel; see [28], [29] and the references therein. These 
works are based on the pure displacement formulation and are concerned mainly with 
compressible materials. 

In this paper, we propose solving the whole indefinite system arising from the mixed 
spectral element discretization using the results in Klawonn [24], [23] and extending his 
h-ve rsion study to spectral elements. We will consider both block-diagonal and block- 
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triangular preconditioners. In the first case, the preconditioned operator is symmetric 
indefinite and we can use the PCR method. In the second case, the preconditioned 
operator is no longer symmetric and we will consider three iterative methods: GMRES 
without restart, Bi-CGSTAB and QMR; see Barret et al. [5] and Freund, Golub, and 
Nachtigal [19] for an introduction to these methods. 

The main result of this paper is that the convergence rate of the proposed algorithms 
is independent of the penalty parameter v, the number of spectral elements N and 
mildly dependent on the the spectral degree n.via the inf-sup constant. This is due 
to the dependence on n of the inf-sup constant for our choices of spectral element 
spaces in the discretization. We will consider two choices of mixed spectral spaces, 
known as the Q n - Q n -i and Q n - P„_i methods; see Maday, Patera and Rpnquist 
[26] and Stenberg and Suri [36]. Several numerical experiments reported in the final 
section confirm this result and show that the number of iterations required by the 
triangular preconditioner is much smaller than the number of iterations required by 
the block-diagonal preconditioner, while its cost is only marginally higher. 

The organization of the paper is as follows. Section 2 introduces the elasticity and 
Stokes systems in both the pure displacement and mixed formulation. In Section 3, two 
mixed spectral element discretizations are introduced. These are based on the Gauss- 
Legendre-Lobatto (GLL) quadrature, briefly reviewed in Section 3.1. The known results 
for the associated inf-sup constants are reported in Section 3.2. The preconditioned 
iterative methods and the main convergence results are introduced in Section 4, with 
the block-diagonal preconditioner in 4.1 and the triangular preconditioner in 4.2. In 
Section 5, we report the results of several numerical experiments in three dimensions, 
both with block-diagonal and triangular preconditioners, with one and many spectral 
elements. 


2 The linear elasticity and Stokes systems 


We consider a polyhedral domain ft C -ft 3 , fixed along a subset of its boundary To, 
subject to a surface force of density g along Ti = dQ — To and subject to an external 
force f. Let V be the Sobolev space V = {v£ if 1 (fl) 3 : v|r 0 = 0}. The linear elasticity 
problem (pure displacement model) consists in finding the displacement u € V of the 
domain $7 such that: 

2 fx I c(u) : c(v) dx + A / div u divv dx = < F, v > Vv 6 V, (1) 

J q Jn 

where A and /z are the Lame constants, c^(u) = + §^-) is the linearized stress 

tensor, and the inner products are defined as 


c(u) : €(v) = 


3 3 


EE f >j( u ) £ ij( v )- 

»= 1 j=l 


< F, v > = 



fiv, dx + 



g t v t ds. 


Almost incompressible materials are characterized by very large values of A, or, in 
terms of the Poisson ratio v = 2(A ^ - , -, by v close to 1/2. When low order h - version 
finite elements are used in the discretization of (1), the locking phenomenon causes a 
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deterioration of the convergence rate as h -* 0; see Babuska and Suri [3]. If the p - version 
is used instead, locking in u is eliminated, but it could still be present in quantities of 
interest such as Xdivu. Moreover, the stiffness matrix obtained by discretizing the pure 
displacement model (1) has a condition number that goes to infinity when v — *• 1/2. 
Therefore, the convergence rate of iterative methods deteriorates rapidly as the material 
becomes almost incompressible. 

These locking problems can be overcome by introducing the new variable p = 
-Xdivu G L 2 (JI) = W and by rewriting the pure displacement problem in the following 
mixed formulation (see Brezzi and Fortin [12]). Find (u,p) G V x W such that 


f 2p f Q e(u) : e(v) dx - j Q div\pdx = <F,v> Vv G V , 

\ ~ So divu q dx - j f Q pq dx = 0 Vq G W. ' 

With the standard definitions a £ (u, v) = 2p f n e(u) : e(v)dx , b(v,q ) = - f Q divvqdx, 
and c(p, q) = f Q pqdx , this problem takes the form: find (u T p) G V x W such that 


a c (u,v) + 6(v,p) = <F,v> Vv G V 

b(u,q) - jc(p,q ) = 0 Vq 6 W. 


(3) 


When A — t oo (or, equivalently, v — * 1/2), we obtain from (2) the limiting problem for 
incompressible elasticity: 


a £ (u,v) + b(v,p) = <F,v> Vv G V 

b(u,q) = 0 Vq£W. 


(4) 


In case of homogeneous Dirichlet boundary conditions on the whole boundary dfi, 
the pressure will have zero mean value, so we define W — Tq(^)- case > problem 

(2) can equivalently be written in the following way (see Brezzi and Fortin [12]): 


g(u,v) + b(v,p) = <F,v> Vv G V 

b(u,q) - j^c(p,q) - 0 Vq G W, 


( 5 ) 


where here a(u, v) = p f n Vu : Vv dx. The limiting problem when A -+ oo is the Stokes 
system describing the velocity u and pressure p of a fluid of viscosity p: 

f a(u, v) + b{v,p ) = < F,v > Vv G V 

\ b(u,q) = 0 Vq€W. KO> 

The penalty term in (5) might be present due to stabilization techniques. 


3 Mixed spectral element methods 

For an introduction to spectral elements see Patera [30], Maday and Patera [27], Maday, 
Patera and Ronquist [26] and the references therein. 

Let Sl re f be the reference cube [-1, l] 3 and let Q n {^ref) be the set of polynomials on 
Sl Te f of degree n in each variable and P n (Q Te j) be the set of polynomials on Q Te j of total 
degree n. Let the domain be decomposed into a finite element triangulation (jflj 
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of nonoverlapping elements. Each ft, is the affine image of the reference cube ft, = 
F,(ft re/ ), where F, is an affine mapping. We discretize each displacement component 
by conforming spectral elements, i.e. by continuous, piecewise polynomials of degree 
n: 

V n = {v € V : Vfcln, o Fi € Q„(ft re /), * = 1, • • - , N, k = 1,2,3}. 

We consider two choices for the discrete pressure space W n : 

W? = {q 6 W : q, o F, 6 g„_ 2 (ft re /), i = 1, ■ • - , JV}, 

W 2 n = {q € W : « o F t 6 P«_i (««/), t = 1, • • • , AT}. 

The choice WJ* has been proposed for the Stokes system by Maday, Patera and Ronquist 
[26] and it is known as the Q n - Q n - 2 method. A basis for W* can be constructed 
by using the tensor-product Lagrangian interpolants associated with the internal GLL 
nodes, described in the next section in more details. 

The second choice corresponds to Method 2 analyzed in Stenberg and Suri [36]. 
We will call this method Q n - P n -i- For P n _ 1 it is not possible to have a tensorial 
basis, but other standard bases, common in the p-version finite element literature, can 
be used. 

3.1 Gauss-Lobatto-Legendre (GLL) quadrature and the discrete prob- 
lem 

The efficient evaluation of the multiple integrals of polynomials, involved in our model 
problem, is based on numerical quadrature at the GLL points. Let {&, 0 

the set of GLL points on the reference cube [-1, l] 3 , and let <r, be the weight associated 
with & . Let l t {x) be the Lagrange interpolating polynomial vanishing at all the GLL 
nodes except at &, where it equals one. By tensor product; the basis functions on the 
reference cube are then defined by 

li(x)lj(y)lk{z), 0 <i,j,k<n. 

Since every polynomial in Q n (Q re j) can be written as 

u(x,y,z) = 

,=0.7=0 k = 0 

these basis functions form a nodal basis. Each integral of the continuous model (2) is 
then replaced by GLL quadrature sums. On ft r€ / 

n n n 

{u,v)Q,n rtJ = 

2—0 j = 0 /c=0 


and in general on fi 

N n 

} ] (ti 0 F s ){v 0 £,k)®i&jGkt 

5=1 tj,k — 0 
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Figure 1: Sparsity pattern of the stiffness matrix K for model problem (7) discretized 
with method Q n - Q n -2 on one element, n = 5, v - 0.3 





where | J s | is the determinant of the Jacobian of F s . The analysis of this discretization 
technique can be found in Bernardi and Maday [7] and Mad ay. Patera and Ronquist 
[26]. 

The discrete problem obtained from pb. (3) is: 


j a e Q ( u,v) + b Q (v,p) = < F,v >Q t n Vv e V" 

\ Mu,?) - \c Q (p,q) = 0 VqeVF", 


where Oq(u,v) = 2//(e(u) : c(v))q, q, M v »? ) = -( divv,q)Q c(p,q ) = (p,?)q, n- 
The bilinear forms 6(-, •) and c(-, •) are computed exactly by. GLL quadrature since the 
ft, are affine images of the reference cube. This system is a saddle point problem with 
a penalty term and has the following matrix form: 


K x = 


A B t 
B -\C 


x = b . 


( 8 ) 


The stiffness matrix K is symmetric and indefinite. It is less sparse than the one 
obtained by low-order finite elements, but is still well-structured. See Figure 1 for the 
sparsity structure of K . An analogous discrete problem with C = 0 is obtained in 
the incompressible case. For the Stokes problem, the discretization of the equivalent 
formulations (5) and (6) lead to an analogous block structure, with A consisting of three 
uncoupled discrete laplacians and with the penalty term in (5) scaled by 1/(A + p). 
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3.2 Estimates of the inf-sup constant for spectral elements 


The convergence of mixed methods depends not only on the approximation properties 
of the discrete spaces V n and but also on a stability condition known as the 

inf-sup (or LBB) condition; see Brezzi and Fortin [12]. For numerical studies of the 
inf-sup constant of various /i- version finite elements, see Bathe and Chapelle [6] and 
Aristov and Chizhonkov [1]. While many important h-version finite elements for Stokes 
problems satisfy the inf-sup condition with a constant independent of h , the important 
spectral elements proposed for Stokes problems, such as the Q n — Q n ~ 2 and Q n — P n ~i 
methods, have an inf-sup constant that approaches zero as ( d = 2,3). This 

result has been proven for the Q n - Q n -2 method by Maday, Patera and Rtfnquist [26], 
where an example is constructed showing that this estimate is sharp. Stenberg and 
Suri [36] proved the following, more general, result covering both methods. 

Theorem 1 ( Stenberg and Suri [36]) Let the spaces V n and W n satisfy assumptions 
(A1)-(A4) of [36] (satisfied by both our methods). Then for d = 2,3 


sup 

veV"\{o} 


(diw, q ) 

iMi^r 


> Cn-^\\q\\ L , 


v<? e w n . 


where the constant C is independent of n and q. 


For the Q n - P n -\ method, no example is known regarding the sharpness of this 
estimate. We can rewrite the inf-sup condition in matrix form as 

q'BA-'B'qyplq'Cq V 9 G W n , (9) 


where So = Cn~( V - ) is the inf-sup constant of the method;, see Brezzi and Fortin [12]. 
Therefore Si scales as A mtn (C _1 BA~ 1 B t ) . If S\ is the continuity constant of B, we 
have 


v ( B ( q < /3 1 ( 9 ‘C' 9 ) 1/2 (v t Av) 1/2 Vv € V n ,V<? € W n . 


( 10 ) 


From (9) and (10) it follows that the 

q t BA~ l B t q 


So < 


q'Cq 


<s\ 


V<? € W n , 


for positive constant So and Si- For stable h- version finite elements, both So and Si 
are independent of h. Theorem 1 shows that this is no longer the case when spectral 
elements are used. However, numerical experiments by Maday, Meiron, Patera and 
Ronquist [25] and [26], have shown that for the Q n - Q n -2 method, for practical values 
of n (e.g. n < 16), the dependence of So on n is much weaker. In our numerical 
experiments in Section 5, we show that the situation is even better for the Q n - P n ~i 
method. Of course, the trade-off in this case is the loss of a tensorial basis. 


4 Preconditioned iterative methods 

The indefinite system Kx = b obtained from our spectral element discretization (8), 
will be solved iteratively by preconditioned Krylov methods for indefinite systems. Two 
classes of preconditioners will be considered: block-diagonal and triangular. 


6 



4.1 Block-diagonal preconditioners 

We first consider a block- diagonal preconditioner for K with positive definite blocks A 
and C: 

D = 

We will denote by D the case with exact blocks A = A and C = C. Interesting choices 
for A are given by h-version finite element discretizations on the GLL mesh or by 
substructuring domain decomposition methods, where ao and a\ have a polylogarithmic 
dependence on the spectral degree n (for the scalar case, see Pavarino and Widlund 
[32] and Casarin [14]). Since the resulting preconditioned system is symmetric, we can 
use the Preconditioned Conjugate Residual Method (PCR); see Ashby, Manteuffel and 
Saylor [2] and Hackbusch [21]. See Elman [16] for a short description of the ORTHMIN 
version of PCR for symmetric indefinite systems. 

In his thesis [24], Klawonn considered low-order finite elements and proved an esti- 
mate for the condition number of D~ X K , under the following assumptions that A and 
C are good preconditioners for A and C respectively: 
i)3co,oi > 0 such that 

GoV ( Av < v'Av < ajv^v Vv 6 V n ; 

n) C is spectrally equivalent to the pressure mass matrix C : 3co,Ci > 0 such that 
c&Cq < q t Cq < ctfCq Vq € W n . 

Theorem 2 (Klawonn [24], pp- 46-47) 

cond(D~ l K) < maa: j a i’ c i} cond{D~ l K) 
min{a^CQ} 

and 

cond{D~ l K) < 

where /3 q is the inf-sup constant of the method and (i\ is the continuity constant of B. 

Clearly, this abstract result can also be applied to high-order elements. Combining 
Theorems 1 and 2, we obtain convergence estimates for both methods we have proposed. 

Theorem 3 If K is the stiffness matrix of the discrete system (7) obtained with either 
the Q n - Q n - 2 or the Q n - P n -\ method and D is the block-diagonal preconditioner 
(11), then 

cond(I)~ 1 Ii) < Cfo 1 = Cn^\ d = 2,3. 

We remark that the number of iterations of the PCR algorithm applied to an indefinite 
system is bounded by the condition number of the system (see Hackbusch [21]). This is 
different from the bounds for conjugate gradient algorithms, where the number of iter- 
ations is bounded by the square root of the condition number of the system. Therefore, 
the number of iterations of our preconditioned algorithm is bounded by Cn^~K 


1/2 + yjfi + 1/4 
-1/2 + + 1/4 



t 



4.2 Triangular preconditioners 


An alternative way to precondition the saddle-point problem (8) is provided by the 
lower and upper triangular preconditioners 


T L 


A 0 


A B t 

B C 

, Tv = 

° C 


( 12 ) 


where A and C are positive definite matrices. Again, we will denote by Tl and Tv 
the case with exact blocks A = A and C = C. Since the resulting preconditioned 
system is no longer symmetric or positive definite, we need to use Krylov methods for 
general nonsymmetric systems. In particular, we will consider three relatively recent 
methods: GMRES, Bi-CGSTAB and QMR; see Barret at al. [5] and Freund, Golub 
and Nachtigal [19]. We remark that each application of the inverse of the triangular 
preconditioners Tl or Tv is only marginally more expensive than the block-diagonal 
preconditioner. In fact, both preconditioners require the solution of a system for A and 
one for C . In addition, the triangular preconditioner requires only one application of 
B (or B *): 


' A 

0 ' 

-1 

u 


A~ l 

0 


u 


A~ l u 

B 

C 


P 


-C~ l B A" 1 

c - 1 


P 


C-'i-BA-'u + p) _ 


In Klawonn [24], it is first proved a bound for the spectrum o of the preconditioned 
operator with exact blocks. The surprising result is that such spectrum is a subset of 
the positive real axis. 

Theorem 4 (Klawonn [24], p. 56) 

o(T[ l K) = <r{KTu l ) C [$,/?? + 1] U {1}. 


Combining this result and the estimate of (3 0 given in Theorem 1 for our spectral 
element spaces, we obtain the following result. 

Theorem 5 If K is the stiffness matrix of the discrete system (7) obtained with either 
Q n _ Q n _ 2 or Q n - P n _ i spectral elements and T is the lower or upper triangular 
preconditioner (12) with exact blocks , then 

cond(T - 1 K) < C(3q 2 = Cn (d ~ l) , d = 2, 3. 


If GMRES is used to solve our problem, it is possible to prove that the number of 
iterations required is bounded by the square root of the condition number of T~ r K; 
see Klawonn [24], Theorems 5.3 and 5.4. Therefore, we have the same bound Cn^~^ 
as for PCR with block-diagonal preconditioner. 

The case of a triangular preconditioner with inexact blocks is studied in Theorem 
5.2 in Klawonn [24], pg. 59, under the standard assumptions i) and ii) of the previous 
section. The estimate provided is analog to the case with exact blocks, but it is more 
complicated and we refer to [24] for the details. 
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5 Numerical results 


The numerical results are divided into a preliminary section regarding the inf-sup con- 
stant and into four sections corresponding to block- diagonal and triangular precondi- 
tioners, each divided into single-element and multi-element case. The iterative meth- 
ods considered are PCR for the block-diagonal preconditioner and GMRES (without 
restart), Bi-CGSTAB and QMR for the triangular preconditioner. All the computations 
were performed in MATLAB 4.2 on a SPARCcenter 2000. The model problem consid- 
ered is (2), discretized with the Q n -Qn- 1 or Q n ~ Pn-i spectral element methods. The 
resulting discrete systems have a structure as in (8). The implementations of GMRES, 
Bi-CGSTAB and QMR are the Matlab templates from [5], while the implementation 
of PCR is the same as in [16]. Except Table 11 showing the discretization errors in the 
i 2 -norm, all the results reported are iteration counts for the iterative methods consid- 
ered. The initial guess is always zero and the right-hand side / consists of uniformly 
distributed random numbers in [-1,1]. The stopping criterion is IIt-,- || 2 / 1| r 0 H 2 < 10 -6 , 
where r; is the i-th residual. We did not try to optimize any of these routines and in 
each table, the size of the largest problem we were able to run was determined by the 
size of the available memory. This was particularly limiting in the multi-element case, 
where already with four elements, we could run only cases up to n — 6. We considered 
only preconditioners with exact blocks, in order to study the algorithms under the best 
of circumstances. For the single-element block- diagonal case, we considered also pre- 
conditioners with inexact blocks based on piecewise linear finite elements on the GLL 
mesh. 

5.1 The inf-sup constant 

We first report in Table 1 a comparison of the spectrum of the matrices C -1 5 = 
C -1 B A~ x B l associated with the two methods Q n ~Qn- 2 and Q n — P n - 1 - Since the inf- 
sup constant /? 0 scales like x/Amtn , these results give an indication on the performance 
of the PCR method reported in the following tables. The first set of results for the 
Qn - Qn- 2 method agree with the 3-d results of Maday, Patera and Rpnquist [26]. For 
these relatively low values of n, A mtn scales like n -1 and therefore /3 0 scales like n -0 - 5 , 
which is better than the value predicted by the theory (n -1 ). 2-d numerical results 
in Maday, Meiron, Patera and Rpnquist [25] for higher values of n (16 < n < 36) 
show that the decay of /5 0 approaches the theoretical bound, but is still better than 
the value predicted by the theory (rc -0 ' 5 ). The case n = 10 could not be run due to 
memory limitations. The second set of results in Table 1 show that the Q n — P n ~\ 
method has a much better inf-sup constant. From so few values of n, it might look 
like A m i n is bounded away from zero. However, a closer look shows that of A, nrrl has 
now a zig-zag behavior. By separating odd and even values of n, we found that A min 
scales approximately like n~ 01 . Higher values of n are needed in order to understand 
the asymptotic behavior of A mm for this method. The maximum eigenvalue quickly 
approaches the same value A maI = 0.65 for both methods. 

We remark that for h-ve rsion finite elements, a numerical study of the inf-sup 
constant seems simpler. In Bathe and Chapelle [6], only three or four values of h are 
needed to predict the asymptotic behavior of A mtn . 
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Table 1: Inf-sup constant (3p - A t7 ; m 


n 

Qn 

cond{C~ l S) 

- Qn-2 

^max 

^min 

Qn 

cond{C~ l S) 

-Pn-l 

^ max 

^mtn 

3 

2.1084 

0.2284 

0.1083 

2.1527 

0.3611 

0.1677 

4 

7.3040 

0.6334 

0.0867 

2.6771 

0.4570 

0.1707 

5 

9.2670 

0.6447 

0.0695 

3.4258 

0.5973 

0.1743 

6 

11.2829 

0.6500 

0.0576 

3.7291 

0.6097 

0.1635 

7 

13.4537 

0.6500 

0.0483 

3.9161 

0.6499 

0.1659 

8 

15.8016 

0.6500 

0.0411 

4.0713 

0.6499 

0.1596 

9 

18.3445 

0.6500 

0.0354 

4.0446 

0.6500 

0.1607 

10 

- 

- 

- 

4.1653 

0.6500 

0.1560 


Table 2: Condition numbers for v = 0.5; exact block- diagonal preconditioner 


n 

Qn-Qn 

cond(D~ l K) 

-2 

cond(K) 

Qn ~ Pn 
cond(D~ 1 K) 

-l ! 

cond(K) 

3 

10.827 

29.338 

7.931 

69.613 

4 

16.259 

130.52 

8.196 

214.63 

5 

20.040 

404.29 

8.556 

869.92 

6 

23.977 

1,042.1 

9.096 

3,656.6 

7 

28.332 

2,371.1 

9.121 

18,427 

8 

33.038 

4,609.2 

9.437 

86,673 

9 

38.132 

8,845.7 

. 9.382 

452,767 
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Table 2 reports the condition numbers of the preconditioned system here 

K is the Stokes matrix obtained for v — 1/2 and D is the preconditioner with exact 
blocks. By Theorem 3, these condition numbers scale like the inverse of the respective 
inf-sup constant. In fact, the results for Q n - Q n - 2 clearly show a linear growth with 
n. The results for Q n - P „_ 1 are much better and, in comparison, they almost look 
bounded by a constant. However, by again separating odd and even values of n, the 
growth still appears linear. 

5.2 Block-diagonal preconditioner: single-element case 

In Table 3, we report the PCR iteration counts for both methods with exact precon- 
ditioner. We followed the PCR implementation of Elman [16], which switches from 
the ORTHOMIN to the ORTHODIR version to avoid breakdown. In our experiments, 
this switch often took place for v near and equal to 1/2. As in Klawonn [22], the 
results are uniform in the Poisson ratio v\ for each fixed degree n, the number of PCR 
iterations is bounded by a constant independent of v. As the material becomes almost 
incompressible, the number of iterations tends to a constant which is the number of 
iterations required by the limiting Stokes problem. As the spectral degree n increases, 
the number of iterations increases, in agreement with Theorem 3. This effect is less 
pronounced for compressible materials (for v = 0.3 and 0.4 the number of iterations 
stays practically constant), but becomes more important near or at the incompressible 
limit. This is particularly true for the Q n - Q n - 2 method, where the growth of the 
number of iterations for v = 1/2 is clearly linear. The results for Q n - P n _i are better, 
as expected from the better inf-sup constant of this method. In this case, it is even 
hard to read a linear growth from the table, which has large constant blocks. Graphs 
showing the convergence history of both methods for n — 8 and u = 1/2 can be found 
in Figure 2. In Table 4, the same results are reported for the equivalent formulation 
(5) instead of (2). This implies that block A in K now consists of three uncoupled 
discrete laplacians, one for each component of u. The problem is somewhat harder to 
solve and PCR takes more iterations than in each corresponding case of the previous 
table (except n = 3 for Q n - P„_ 1 ). Again, the results for Q n - P„_i are consistently 
better than those for Q n —Q n - 2 - Now a linear growth with n is clear for both methods 
(for Q n - P n _ 1 the odd and even values of n have to be separated). 

Next, we consider an inexact preconditioner by choosing as u-block A the Q\ finite 
element stiffness matrix obtained by discretizing the term f n Vu : Vv dx on the GLL 
grid. In the scalar case, it is well-known that such matrix is spectrally equivalent to 
the stiffness matrix obtained by spectral discretization; see Deville and Mund [15]. In 
Table 5, we study numerically the quality of such preconditioner in three dimensions 
for a Poisson problem on the reference cube with homogeneous Dirichlet boundary con- 
ditions. In the first column, we report the condition number of Fq*K. It is not obvious 
that the values are bounded by a constant, but they appear to grow slower then any 
power of n or log(n), as results from log-log plots. In any case, these values are larger 
than the corresponding ones reported by Rpnquist [33] for Fp, (the Pi finite element 
stiffness matrix obtained by dividing each element of the GLL mesh into tetrahedra). 
The values reported by Ronquist are all from 2 to 2.65, for values of n between 4 and 
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Table 3: Iteration counts with exact block-diagonal preconditioner 






Qn 

— Qn- 2 




n 





V 





0.3 

0.4 

0.49 

0.499 

0.4999 

0.49999 

0.499999 

0.5 

2 

1 

1 

1 

1 

1 

1 

1 

1 

3 

7 

7 

7 

7 

7 

7 

7 

7 

4 

ii 

13 

21 

21 

21 

21 

21 

21 

5 

ii 

15 

27 

31 

31 

31 

31 

31 

6 

ii 

15 

29 

33 

35 

35 

35 

35 

7 

n 

15 

29 

35 

35 

35 

37 

37 

8 

ii 

15 

31 

37 

39 

39 

39 

39 

9 

ii 

15 

31 

39 

41 

41 

41 

41 





Qn 

-Pn-i 




n 





V 





0.3 

0.4 

0.49 

0.499 

0.4999 

0.49999 

0.499999 

0.5 

2 

2 

2 

2 

2 

2 

2 

2 

2 

3 

9 

9 

9 

9 

9 

9 

9 

9 

4 

9 

11 

15 

15 

15 

15 

15 

15 

5 

9 

13 

19 

21 

21 

21 

21 

21 

6 

9 

13 

21 

21 

21 

21 

21 

21 

7 

11 

13 

21 

21 

21 

21 

21 

21 

8 

11 

13 

23 

25 

25 

25 

25 

25 

9 

11 

13 

23 

25 

25 

25 

25 

25 

10 

11 

13 

23 

25 

25 

25 

25 

25 
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Table 4: Equivalent formulation: iteration counts with exact block- diagonal precondi- 
tioner , 






Qn 

“ Qn-2 




n 

0.3 

0.4 

0.49 

0.499 

V 

0.4999 

0.49999 

0.499999 

0.5 

2 

1 

1 

1 

1 

1 

1 

1 

1 

3 

7 

7 

7 

7 

7 

7 

7 

7 

4 

15 

17 

21 

23 

23 

23 

23 

23 

5 

19 

21 

33 

37 

37 

37 

37 

37 

6 

17 

21 

37 

41 

43 

43 

43 

43 

7 

19 

21 

39 

45 

45 

45 

45 

45 

8 

19 

23 

41 

47 

47 

47 

47 

47 

9 

19 

21 

41 

47 

49 

49 

49 

49 





Qn 

- Pn-i 




n 





V 





0.3 

0.4 

0.49 

0.499 

0.4999 

0.49999 

0.499999 

0.5 

2 

2 

2 

2 

2 

2 

2 

2 

2 

3 

7 

7 

7 

7 

7 

7 

7 

7 

4 

10 

13 

15 

15 

15 

15 

15 

15 

5 

11 

15 

23 

25 

25 

25 

25 

25 

6 

13 

17 

25 

27 

27 

27 

27 

27 

7 

13 

17 

27 

29 

29 

29 

29 

29 

8 

13 

17 

31 

33 

33 

33 

33 

33 

9 

13 

17 

29 

31 

31 

31 

31 

31 

10 

13 

18 

31 

33 

35 

35 

35 

35 


Table 5: Q\ finite element preconditioner on the GLL mesh for Poisson equation: 
condition numbers and relative errors with known exact solution 


n 


BBS1 

■HBH 

3 

4.8150 

0.9794 

0.0130 

4 

8.4566 

0.5892 

0.0020 

5 

11.1569 

0.3481 

5.1163e-5 

6 

13.0747 

0.2352 

1.4581e-5 

7 

14.4623 

0.1704 

2.1302e-7 

8 

15.4977 

0.1296 

7.4257e-8 

9 

16.2954 

0.1021 

7.0685e-10 

10 

16.9275 

0.0826 

2.7414e-10 

11 

17.4406 

0.0682 

1.8431e-12 

12 

17.8653 

- 

- 


13 





14 






Figure 2: Convergence history for n = 8 and v — 1/2: QP — ex = Q n - P n _i method 
with exact preconditioner, QQ — ex = Q n — Q n - 2 method with exact preconditioner, 
QP - Ql = Q n - P n - 1 method with inexact Q\ preconditioner, QQ - Ql = Q n - Q n - 2 
method with inexact Q\ preconditioner 



12. In the second and third columns of Table 5, we report the relative errors in the 
discrete / 2 -norm between the exact solution u ex = sin |(x + 1) sin |(y + 1) sin j(z + 1), 
obtained by computing the appropriate right-hand side / = -Au, and the discrete so- 
lution u n (spectral) or uq 1 ( Qi fern on the GLL mesh). The difference between spectral 
and h-ve rsion finite element accuracy is very clear. 

In Table 6, we report the iteration counts for the model problem (2) when the 
inexact block A — diag(FQ 1 , Fq 1 ^ Fq^ is used in the preconditioner. Even if the 
uniformity in v is preserved, the number of iterations grows considerably, especially for 
higher values of n. Therefore, it is does not appear that this inexact preconditioner is 
effective for PCR methods applied to mixed spectral systems. Results for the equivalent 
model problem (5) were similar and are not reported. 

Figure 2 shows the convergence history of the Q n - Q n - 2 and Qn ~ Pn - 1 methods 
with n = 8 and with exact and inexact Q\ preconditioners for the Stokes problem. The 
resulting graphs are similar to the ones reported in Elman [16]. 

5.3 Block-diagonal preconditioner: multi-element case 

Tables 7 and 8 report the iteration counts for Q n -Q n - 2 and Q n -P n - 1 respectively. 
Here we study the dependence of the number of iterations on the number of elements N 
for a fixed spectral degree n. This is analog to studying the dependence on h for a ftp- 
finite element method. We divide the domain Q into jV = N x x N y x N- subcubes and 
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Table 7: Iteration count for Q n ~Qn-2 with exact block-diagonal preconditioner: chang- 
ing N for fixed n 


n 

N = N x x N y x N z 

0.3 

0.4 

0.49 

0.499 

V 

0.4999 

0.49999 

0.499999 

0.5 

2 

1 = 1 x 1 x 1 

1 

1 

1 

1 

1 

1 

1 

1 

2 

8 = 2 x 2 x 2 

6 

7 

7 

7 

7 

7 

7 

7 

2 

27 = 3 x 3 x 3 

9 

11 

17 

17 

17 

17 

17 

17 

2 

64 = 4 x 4 x 4 

10 

13 

19 

21 

21 

21 

21 

21 

2 

125 = 5 x 5 x 5 

10 

13 

21 

23 

23 

23 

23 

23 

3 

1 = 1 x 1 x 1 

7 

7 

7 

7 

7 

7 

7 

7 

3 

8 = 2 x 2 x 2 

11 

13 

21 

23 

23 

23 

23 

23 

3 

27 = 3 x 3 x 3 

11 

13 

21 

23 

23 

23 

23 

23 

3 

64 = 4 x 4 x 4 

11 

13 

21 

23 

25 

25 

25 

25 

4 

1 = 1 x 1 x 1 

11 

13 

21 

21 

21 

21 

21 

21 

4 

8 = 2 x 2 x 2 

11 

15 

23 

27 

27 

27 

27 

27 


Table 8: Iteration count for Q n -P n - \ with exact block-diagonal preconditioner: chang- 
ing N for fixed n 


n 

< 

II 

fei 

x ^ ^ y ^ ^ z 





V 



M 



0.3 

0.4 

0.49 

0.499 

0.4999 

0.49999 

0.499999 

0.5 

2 

1 = 

lxlxl 

2 

2 

2 

2 

2 

2 

2 

2 

2 

8 = 

2x2x2 

9 

13 

19 

21 

21 

21 

21 

21 

2 

27 = 

3x3x3 

11 

14 

23 

27 

27 

27 

27 

27 

2 

64 = 

4x4x4 

11 

15 

23 

27 

27 

27 

27 

27 

2 

125 = 

=5x5x5 

11 

15 

23 

27 

27 

27 

27 

27 

3 

1 = 

lxlxl 

9 

9 

9 

9 

9 

9 

9 

9 

3 

8 = 

2x2x2 

11 

13 

21 

23 

23 

23 

23 

23 

3 

27 = 

3x3x3 

11 

13 

23 

25 

25 

25 

25 

25 

4 

1 = 

lxlxl 

9 

11 

15 

15 

15 

15 

15 

15 

4 

8 = 

2x2x2 

11 

13 

23 

25 

25 

25 

25 

25 
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Table 9: Iteration count for Q n —Qn- 2 with exact block-diagonal preconditioner: chang- 
ing n for fi xed JV = 4 = 2x2xl 


71 

0.3 

0.4 

0.49 

0.499 

V 

0.4999 

0.49999 

0.499999 

0.5 

2 

5 

5 

5 

5 

5 

5 

5 

5 

3 

10 

13 

21 

23 

23 

23 

23 

23 

4 

11 

15 

25 

29 

29 

29 

29 

29 

5 

11 

15 

27 

31 

33 

33 

33 

33 

6 

11 

15 

27 

35 

35 

35 

35 

35 


Table 10: Iteration count for Q n - P n -i with exact block- diagonal preconditioner: 
changing n for fixed N = 4— 2x2x1 


71 

0.3 

0.4 

0.49 

0.499 

V 

0.4999 

0.49999 

0.499999 

0.5 

2 

9 

11 

17 

19 

19 

19 

19 

19 

3 

10 

13 

19 

23 

23 

23 

23 

23 

4 

11 

14 

23 

25 

27 

27 

27 

27 

5 

11 

14 

23 

27 

27 

27 

27 

27 

6 

11 

14 

23 

27 

27 

27 

27 

27 


we take N x = N y = N z in order to always have a cubic domain (and avoid comparing 
problems with different aspect ratios). Due to the cubic growth of N , we could only 
run cases with a low value of n = 2,3,4. The results of the tables seem to indicate a 
bound on the number of iterations independent of TV, in agreement with the theory. 
This is particularly evident for the Q n - P n _i method and for n = 2, which allows us 
to run with sufficiently many elements. As before, the results are uniform in v and the 
incompressible limit is the hardest case, for each n and N fixed. 

In tables 9 and 10, we study the dependence of the number of iterations on n, for 
a small fixed number of elements N = 4, with N x = N y = 2 and N z = 1. The linear 
growth of the number of iterations with n is clearly visible in the incompressible limit 
for Q n — Qn - 2 (Table 9), while for compressible materials ( v « 0.3 — 0.4) the number 
of iterations seems insensitive to n and bounded by a constant. For the Q n - P n -\ 
method, the results of Table 10 are better and the number of iterations seems bounded 
by a constant also in the incompressible limit. However, as was shown in the single- 
element case, higher values of n might reveal a growth which is still linear, just with a 
better constant in front of the linear term. 
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Table 11: X 2 -errors with known exact solution for v = 0.3; Q n - Q n - 2 with exact 


N = N x x N y x N z 

n 

G 

iter. 

B 

Q 

||U-gmre5 |2 

\\Pgrnres-p_exj\2 

l|u«lk 

l|P«l|2 


2 

1 

1 

1 

4.6502e-l 

1.0000e+l 


3 

2 

1 

2 

1.0027e-l 

4.0346e-l 


4 

4 

3 

4 

1.8246e-2 

7.6835e-2 

1 = 1 X lx 1 

5 

6 

4 

6 

3.7748e-3 

2.5745e-2 


6 

7 

4 

7 

3.6803e-4 

3.0630e-3 


7 

7 

4 

7 

5.9270e-5 

6.8210e-4 


8 

7 

4 

7 

4.0743e-6 

5.0991e-5 


9 

7 

4 

7 

5.3355e-7 

8.7179e-6 


2 

2 

1 

2 

9.4308e-2 

2.5707e-l 

8 = 2 x 2x 2 

3 

6 

4 

6 

. 1.2310e-2 

1.1750e-l 


4 

7 

4 

7 

1.1785e-3 

1.3814e-2 


5 

7 

4 

7 

9.1187e-5 

1.6598e-3 


5.4 Triangular preconditioner: single-element case 

In the following tables, we have used the convention G = GMRES, B = Bi-CGSTAB, 
Q = QMR. In all cases, we have used the (left) lower-triangular preconditioner T L with 
exact blocks. 

In the first part of Table 11, we report the errors in the T 2 -norm between the 
Qn - Qn- 2 spectral element solution and the known exact solution U\ - u 2 = u 3 - 
smi^xjsinij^yjsmi^z), p = Xdiv(u) for the elasticity problem with v = 0.3 on 
the reference cube (i.e. N x = N y = N z = 1). On each row, corresponding to each 
value of n, we report the iteration counts for the three methods and the errors for 
the displacement u and for p. The results clearly show the spectral convergence of 
the discrete solution to the exact solution. The second part of the table shows the 
same results for a multi-element case with 8 elements. Tables 12 and 13 report the 
iteration counts for Q n - Q n - 2 and Q n - P n - 1 respectively, on the reference element. 
For each value of n and u, the results for GMRES, Bi-CGSTAB and QMR are reported. 
As in the block-diagonal case, the results are uniform in v , i.e. for each fixed n, the 
number of iterations tends to the number of iterations of the limiting incompressible 
case. Moreover, for each fixed v , the number of iterations grows at worst linearly 
with n, in agreement with the theory. This is clear at the incompressible limit, while 
away from it the results are much better and in practice bounded independently of 
n. Among the three iterative methods, Bi-CGSTAB requires the least number of 
iterations, in some cases half of those required by GMRES, but it requires twice as 
many applications of the matrix and the preconditioner. Moreover, Bi-CGSTAB shows 
a more irregular convergence behavior than the other two methods. QMR has iterations 
counts in between Bi-CGSTAB and GMRES, often closer to the last one. QMR also 
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Table 12: Iteration count for Q n - Q n - 2 on one element, with exact lower-triangular 


preconditioner; G=GMRES, B=Bi-CGSTAB, Q— QMR 


71 

0.3 

G B Q 

0.4 

G B Q 

0.49 
G B Q 

0.499 
G B Q 

V 

0.4999 
G B Q 

0.49999 
G B Q 

0.499999 
G B Q 

0.5 

G B Q 

2 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

3 

4 2 4 

4 3 4 

4 3 4 

4 3 4 

4 .3 4 

4 3 4 

4 3 4 

4 3 4 

4 

6 3 6 

747 

11 7 11 

11 8 11 

11 10 12 

11 10 12 

11 10 12 

11 10 12 

5 

7 4 6 

8 5 8 

15 9 14 

18 17 18 

19 17 18 

19 14 18 

19 18 18 

19 14 18 

6 

7 4 6 

9 5 8 

16 12 15 

20 17 19 

21 15 20 

21 15 20 

21 15 20 

21 15 20 

7 

7 4 6 

9 5 8 

17 14 16 

22 17 20 

22 17 21 

22 16 21 

22 16 21 

22 20 21 

8 

7 4 7 

9 5 8 

20 14 18 

24 20 21 

24 24 22 

24 20 22 

24 20 22 

24 20 22 

9 

7 4 7 

9 5 8 

20 13 18 

25 19 22 

26 20 23 

26 21 23 

26 19 23 

26 17 23 


Table 13: Iteration count for Q n — P n -i on one element, with exact lower-triangular 
preconditioner; G=GMRES, B-Bi-CGSTAB, Q=QMR 


n 

0.3 

G B Q 

0.4 

G B Q 

0.49 
G B Q 

0.499 
G B Q 

V 

0.4999 
G B Q 

0.49999 
G BQ 

0.499999 

GBQ 

0.5 

GBQ 

2 

2 1 2 

2 1 2 

2 1 2 

2 1 2 

2 1 2 

2 1 2 

2 1 2 

2* 2 

3 

5 3 5 

535 

54 5 

54 5 

54 5 

545 

545 

54 5 

4 

535 

6 4 6 

8 6 8 

86 8 

86 8 

8 6 8 

86 8 

86 8 

5 

6 3 6 

7 4 7 

11 6 10 

11 7 11 

12 7 11 

12 7 11 

12 7 11 

12 7 11 

6 

7 3 6 

8 4 7 

12 7 11 

13 8 11 

13 8 11 

13 8 11 

13 8 11 

13 8 11 

7 

7 4 6 

8 4 7 

13 6 11 

14 7 12 

14 7 12 

14 7 12 

14 7 12 

14 7 12 

8 

7 3 6 

9 4 8 

14 7 12 

15 8 13 

15 8 13 

15 8 13 

15 8 13 

15 8 13 

9 

8 4 6 

9 4 8 

14 7 12 

15 8 13 

15 -8 13 

15 8 13 

15 8 13 

15 8 13 

10 

8 4 6 

10 5 8 

15 7 12 

18 8 13 

18 8 13 

18 8 13 

18 8 13 

18 8 13 


requires twice as many matrix and preconditioner applications compared with GMRES. 
Of course, GMRES without restart requires much more memory than the other two 
methods. In comparison with the block- diagonal results, the triangular preconditioner 
requires many less iterations, sometimes half of those required by PCR. 

5.5 Triangular preconditioner: multi-element case 

Tables 14 and 15 are the analog for triangular preconditioners of Tables 7 and 8 for 
block-diagonal preconditioners. Here, we verify numerically the independence of the 
iteration counts on N while keeping n fixed. Again, we could run with many elements 
only for small values of n = 2,3,4, and we used decompositions of the cubic domain 
into cubic powers of subdomains. The results indicate in all cases an upper bound inde- 
pendent of N. Regarding the different convergence performance of the three methods, 
the same considerations as for the single-element case apply. 
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Table 14: Iteration count for Q n - Q n -2 with exact lower-triangular preconditioner: 


changing N_ for fixed n 


n 

O o 

N 

0.3 

G B Q 

0.4 

G B Q 

0.49 
G B Q 

0.499 
G B Q 

V 

0.4999 
G B Q 

0.49999 
G B Q 

0.499999 
G B Q 

0.5 

G B Q 

2 

~w 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

1 1 1 

2 

2 3 

3 2 4 

4 2 4 

43 4 

4 3 4 

4 3 4 

4 3 4 

4 3 4 

4 3 4 

2 

3 3 

5 3 6 

6 4 7 

85 9 

8 6 9 

8 6 9 

86 9 

8 6 9 

8 6 9 

2 

4 3 

5 4 6 

6 4 7 

9 6 10 

10 7 11 

10 7 11 

10 7 11 

10 7 11 

10 7 11 

2 

5 3 

5 4 6 

7 5 7 

10 7 11 

10 7 12 

11 7 12 

11 7 12 

11 7 12 

11 7 12 

3 


4 2 4 

4 34 

4 3 4 

4 3 4 

4 3 4 

4 3 4 

4 3 4 

4 3 4 

3 

2 3 

6 3 6 

7 4 8 

11 7 11 

12 8 12 

12 8 12 

12 8 12 

12 8 12 

12 8 12 

3 

3 3 

6 4 6 

7 5 8 

11 7 12 

12 8 12 

12 8 13 

12 9 13 

12 9 13 

12 9 13 

4 


6 3 6 

7 4 7 

11 7 11 

11 8 11 

11 10 12 

11 10 12 

11 10 12 

11 10 12 

4 

2 3 

6 4 6 

8 5 8 

12 8 13 

14 9 14 

14 10 15 

14 10 15 

14 10 15 

14 10 15 


Tables 16 and 17 are the analog for triangular preconditioners of Tables 9 and 10. 
Here, we fix a small number of elements N - 4 and we study the iteration counts 
by increasing n and v. Again, the results are uniform in v and linear (at worst) 
in n, with the incompressible case being the hardest one. In comparison with the 
block-diagonal preconditioner, the triangular preconditioner considerably decreases the 
number of iterations, sometimes by as much as one half. 


6 Conclusions 

We have proposed and analyzed iterative methods for the sparse indefinite systems 
arising from the mixed spectral element discretization of elasticity and Stokes prob- 
lems. These systems are solved with a preconditioned conjugate residual method when 
a block-diagonal preconditioner is used or with Krylov methods for nonsymmetric sys- 
tems such as GMRES. Bi-CGSTAB and QMR when a triangular preconditioner is used. 
We have proven and have numerically shown that such algorithms have convergence 
rates bounded by the inverse of the inf-sup constant and independent of the penalty 
parameter in the saddle point formulation (the Poisson ratio for elasticity or a stabi- 
lization parameter for Stokes). The two mixed spectral methods considered, Q n ~Qn- 2 
and Q n - P n - i, have equivalent theoretical convergence bounds, but we have numeri- 
cally shown that the latter one has a better inf-sup constant and gives better iteration 
counts. On the other hand, P n - j does not have a tensorial basis. The exact blocks 
in the preconditioners could be replaced by appropriate preconditioners based on low- 
order discretizations on the GLL mesh and/or domain decomposition techniques. The 
inexact preconditioner based on Q i finite elements on the GLL mesh largely increases 
the iteration counts. Future work should address other inexact preconditioners, such 
as Pi finite elements on the GLL mesh, multigrid or domain decomposition methods 
and preconditioners for the mass matrix C . 
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Table 15: Iteration count for Q n — P n -\ with exact lower-triangular preconditioner: 
changing N for fixed n 


n 

N 

0.3 

G B Q 

0.4 

G B Q 

0.49 
G B Q 

0.499 
G B Q 

V 

0.4999 
G B Q 

0.49999 
G B Q 

0.499999 
G B Q 

0.5 

G B Q 

2 

i 3 

2 1 2 

2 1 2 

■SSI 

2 1 2 

2 1 2 

2 1 2 

2 1 2 

2*2 

2 

2 3 

5 3 6 

6 4 7 


10 7 11 

10 7 11 

10 7 11 

10 7 11 

10 7 11 

2 

3 3 

64 6 

7 5 8 

11 8 12 

13 9 14 

13 9 14 

13 9 14 

13 9 14 

13 9 14 

2 

4 3 

6 4 6 

7 5 8 

11 8 12 

13 9 14 

13 9 14 

13 9 14 

13 9 14 

13 9 14 

2 

5 3 

6 4 6 

7 5 8 

12 8 13 

13 9 14 

13 10 14 

13 10 14 

13 10 14 

13 10 14 

3 

l 3 

53 5 

5 3 5 

5 4 5 

5 4 5 

54 5 

5 4 5 

5 45 

5 4 5 

3 

2 3 

6 4 6 

7 5 7 

10 7 11 

11 7 12 

11 7 12 

11 7 12 

11 7 12 

11 7 12 

3 

3 3 

6 4 6 

7 5 8 

11 7 12 

12 8 13 

12 8 13 

12 8 13 

12 8 13 

12 8 13 

4 

l 3 

5 3 5 

6 4 6 

8 6 8 

8 6 8 

86 8 

86 8 

8 6 8 

8 6 8 

4 

2 3 

6 4 6 

8 5 8 

12 7 12 

12 8 13 

12 8 13 

12 8 13 

12 8 13 

12 8 13 


Table 16: Iteration count for Q n — Q n -2 with exact lower- triangular preconditioner: 
c hanging n for fixed jV = 4 = 2x2xl 


n 

0.3 

G B Q 

0.4 

G B Q 

0.49 
G B Q 

0.499 
G B Q 

V 

0.4999 
G B Q 

0.49999 
G B Q 

0.499999 
G B Q 

0.5 

G B Q 

2 

3 2 3 

3 2 3 

ffFfl 

3 2 3 

3 2 3 

3 2 3 

3 2 3 

3 2 3 

3 

6 3 6 

6 4 7 

BBS 

11 8 12 

11 8 12 

12 812 

12 8 12 

12 8 12 

4 

6 4 6 

8 5 8 

13 9 13 

15 11 15 

15 11 16 

15 13 16 

15 13 16 

15 13 16 

5 

7 4 6 

8 5 8 

15 9 14 

18 14 18 

19 14 19 

19 14 19 

19 14 19 

19 14 19 

6 

7 4 6 

8 5 8 

16 10 15 

20 14 19 

21 16 20 

21 15 20 

21 15 20 

21 15 20 


Table 17: Iteration count for Q n — P n -\ with exact lower-triangular preconditioner: 


ch anging n for fixed A r = 4 = 2x2xl 



0.3 

G B Q 

0.4 

G B Q 

0.49 
G B Q 

0.499 
G B Q 

V 

0.4999 
G B Q 

0.49999 
G B Q 

0.499999 
G B Q 

0.5 

G B Q 

2 

53 5 

6 46 

8 6 9 


9 7 10 



9 7 10 

3 

64 6 

7 4 7 


11 8 12 

11 8 12 

11 8 12 

11 8 12 

11 8 12 

4 

6 4 6 

8 5 8 

12 8 12 

14 8 14 

14 8 14 

14 9 14 

14 8 14 

14 9 14 

5 

64 6 

8 5 8 

13 9 13 

15 12 14 

15 14 14 

15 12 14 

15 13 14 

15 14 14 

6 

7 4 6 

8 5 8 

14 9 13 

15 11 14 

15 11 14 

15 11 14 

15 11 14 

15 11 14 
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